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I.  Background 


This  investigation  is  motivated  by  developing  physical  understanding  of  fundamental  aspects  of 
the  aerodynamics  of  Unmanned  Micro  Air  Vehicles  (UMAVs).  For  a  typical  UMAV  (Wilson, 
1998),  the  small  size  of  the  device  (wing  chord  of  order  2  cm)  and  relatively  slow  flying  speed 
(of  order  10  m/s)  result  in  a  wing  chord  Reynolds  number  Rec  of  order  104.  Understanding  of  the 
aerodynamics  at  this  low  range  of  Reynolds  numbers  is  lacking  as  traditional  aerodynamics 
knowledge  is  based  on  high  Reynolds  number  airfoils  with  thin  attached  boundary  layers  that 
operate  under  steady  or  quasi-steady  conditions.  Faced  by  the  challenges  of  adequate  lift 
generation  in  low  Reynolds  number  flight  based  on  traditional  aerodynamics,  and  observation  of 
nature’s  flying  animals,  a  major  departure  in  UMAV  aerodynamics  is  the  considerations  of 
flapping  wing  designs  where  highly  unsteady  operating  conditions  are  to  be  exploited  instead  of 
avoided  (Lai  and  Platzer,  1999,  Young  and  Lai,  2004,  and  Anderson  et  al.,  1998) .  Added  to  this 
are  the  additional  complexities  introduced  by  highly  3-D  low-aspect  ratio  wings  that  might  also 
include  surface  flexibility,  along  with  high-amplitude  complicated  motion  trajectories,  and  gust 
response/tolerance.  The  problem  is,  indeed,  multi-dimensional  and  immensely  challenging. 

The  behavior  of  flow  around  airfoils  executing  highly  unsteady  motions  has  of  course  been 
studied  for  decades.  The  classical  unsteady  aerodynamic  theory  of  oscillating  airfoils 
(Theodorsen,  1935,  and  von  Karman  and  Sears,  1938)  was  developed  as  a  result  of  interest  in 
aircraft  flutter  problems.  This  theory  later  found  extensive  use  in  bio-fluid-dynamics,  and  the 
propulsion  of  certain  species  of  birds,  insects,  and  aquatic  animals  (Lighthill,  1975,  and  Wu, 
1971).  The  main  elements  of  the  classical  analysis  are  the  potential  flow  along  with  linearized 
boundary  conditions,  small  perturbation  velocities,  and  the  assumption  of  a  planar  vortex  wake. 
All  of  these  assumptions,  appropriate  for  thin  airfoils  under  small  amplitude  unsteady  motions 
along  with  thin  boundary  layers  over  the  airfoil,  fail  for  the  low  Reynolds  number  unsteady 
airfoils  of  UMAVs.  Even  though  nonlinear  processes  such  as  rolled-up  wake  vorticity  have  been 
addressed  using  numerical  techniques  (Giesing,  1968,  Djojodihardjo  and  Widnall,  1969,  and 
Vezza  and  Galbraith,  1985),  all  of  these  approaches  still  rely  on  potential  flow  theory  and  vortex 
methods  where  the  strength  of  the  vorticity  shed  into  the  wake  is  determined  from  the  application 
of  the  Kutta  condition  at  the  trailing  edge.  In  low  Reynolds  number  airfoils  where  the  boundary 
layer  thickness  at  the  trailing  edge  can  be  significantly  larger  than  the  trailing  edge  thickness,  and 
under  highly  unsteady  conditions,  the  conditions  that  should  be  imposed  at  the  trailing  edge 
become  less  clear  than  the  Kutta  condition  used  in  inviscid  analysis  (Koochesfahani  and  Bohl, 
2002). 

Relevant  to  current  UMAV  efforts  is  also  the  large  body  of  knowledge  that  is  available  from 
previous  studies  connected  to  dynamic  stall  (McCroskey,  1982)  and  significant  developments  in 
computational  capabilities  for  such  high  amplitude  unsteady  airfoil  motions,  in  particular  those  at 
AFRL  which  are  of  interest  to  this  proposal  (Visbal,  1986,  and  Visbal  and  Shang,  1989). 
Nevertheless,  recent  studies  in  UMAV  aerodynamics  have  lead  to  the  conclusion  that  the  current 
state  of  knowledge,  predictive  capabilities,  and  experimental  data  are  insufficient  regarding  the 
fundamental  unsteady  aerodynamics  of  low  Reynolds  number  UMAVs  (Shyy  et  al.,  2008) .  The 
coordinated  experimental  and  computational  efforts  undertaken  in  the  present  investigation  is 
aimed  at  improving  our  fundamental  understanding  of  low  Reynolds  number  unsteady 
aerodynamics  and,  in  conjunction  with  other  AFOSR  supported  efforts  and  AFRL,  contributing 
to  improved  predictive  capabilities. 
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II.  Objectives 


The  specific  focus  of  the  present  project  is  the  unsteady  aerodynamics  of  oscillating  airfoils 
with/without  a  flexible  trailing  edge  in  the  chord  Reynolds  number  range  of  2xl03  -  6xl04.  The 
study  encompasses  various  airfoil  motion  trajectories  (pure  pitch,  pure  plunge,  combined  pitch- 
plunge,  asymmetric  trajectory,  etc.),  and  a  range  of  oscillation  frequencies  and  amplitudes.  For 
such  airfoils,  we  seek  to  establish  the  connection  between  the  airfoil  motion  trajectory,  the 
trailing  edge  flexure,  the  time  history  of  vorticity  flux  at  the  trailing  edge,  the  pattern  of  shed 
vorticity  and  its  evolution,  and  the  conditions  for  generation  of  thrust/drag  and  lift.  In  addition, 
we  examine  how  these  processes  are  influenced  by  aeroelastic  parameters  of  the  trailing  edge 
and  aspect  ratio.  The  investigation  is  carried  out  using  coordinated  experiments  and 
computations. 

III.  Methods 

IV.  1.  Three-Degree  of  Freedom  (3  DOF)  System  and  Water  Tunnel  Modification 

The  proposed  work  required  significant  infrastructure  modification  of  the  61  cm  x  61  cm  water 
tunnel  test  facility  at  the  Turbulent  Mixing  and  Unsteady  Aerodynamics  Laboratory  (TMUAL)  at 
MSU.  This  modification,  which  was  completed  successfully,  included  two  parts:  (I)  the  addition 
of  two  large  quartz  inserts  (45.7  cm-high  x  89  cm-wide  x  12.7  cm-deep)  into  one  of  the  acrylic 
side  walls  of  the  61  cm  x  61  cm  x  243  cm  test  section  (to  allow  UV  laser  light  to  pass  through 
for  MTV  measurements);  (II)  the  acquisition,  assembly,  mounting  and  operation  verification  of  a 
three  degree  of  freedom  (3  DOF)  rig  to  allow  pitch,  heave,  and  perch  maneuvers  of  rigid  and 
flexible  airfoils.  Details  concerning  these  modifications  may  be  found  in  Figures  1  and  2  and 
associated  captions. 


3  DOF  motion  system 


Test  section 


Quartz  panels 


UV  mirror 


UV  laser  beam 


Figure  1.  Water  tunnel  test-section  cut-out  showing  the  new  quartz  inserts  and  the  3  DOF  motion  system. 
The  sketch  also  depicts  a  typical  UV  laser  beam  path  for  Molecular  Tagging  Velocimetry  measurements 
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Linear  motor  (spanwise  -  plunge): 
Parker  Trilogy  T4D 


Linear  motor 
(streamwise) 


Six-component  balance: 
ATI  SI-145-5 


Rotary  motor  (pitch): 
Parker  MPP 


Figure  2.  Details  of  the  3  DOF  Motion  System 


IV.  2.  Airfoil  Fabrication  Using  3D  Printing  Methods 

For  the  parts  of  the  investigation  focused  on  studying  the  effect  of  trailing  edge  flexibility, 
modem  multi-material  3D  printing  technology  is  employed  to  fabricate  the  airfoil  models.  With 
this  approach,  the  entire  airfoil  is  fabricated  using  a  3D  printer  (Connex350)  which  is  capable  of 
printing  objects  using  two  materials,  one  rigid  (YeroWhitePlus)  and  one  flexible  (TangoBlack), 
resulting  in  seamless  integration  of  the  rigid  and  flexible  parts  of  the  airfoil.  Two  different  airfoil 
shapes  are  considered:  a  “long-tailed”  and  NACA-0009  airfoil  (see  Figure  3).  The  long-tailed 
airfoil  consists  of  a  rigid  “head”  portion  that  has  a  NACA  0036  profile,  and  a  “tail”  piece  that  is 
made  from  either  the  rigid  or  the  flexible  material.  The  chord  of  the  airfoil  is  c  =  30  mm,  with  the 
head  stretching  over  7.5  mm  length.  The  airfoil’s  thickness  ratio  based  on  c  is  0.09.  To  vary  the 
flexibility  of  the  airfoil,  the  thickness  t  of  the  tail  is  varied.  The  non-dimensional  parameter  G\ 
typically  used  to  characterize  the  structure  flexibility  is  defined  as: 

^  _Et3/(  12(1 -v2)) 

- 

pK4  [l] 

where  E  and  vare  Young’s  modulus  and  Poisson’s  ratio  for  the  tail  material  respectively,  p is  the 
fluid  density,  LA  is  the  ffeestream  velocity  and  C2  is  the  length  of  the  tail  (22.5  mm). 
Measurements  conducted  with  this  airfoil  are  done  in  %  scale  water  tunnel  having  a  test  section 
with  a  flow  area  of  15.2  cm  x  15.2  cm.  With  the  airfoil  spanning  the  width  of  the  test  section,  the 
airfoil’s  aspect  ratio  is  approximately  5. 


3 


Head:  NACA  0036  -  rigid 


Tail:  variable  flexibility 


■  I 

■*-C)=c/4— ■ «j« - c2 -  »j 

* - c  =  30  mm - - 


Figure  3.  Illustration  of  the  airfoil  shapes  employed  in  the  unsteady  airfoil  studies  in  the  present  work: 
long-tailed  (left)  and  NACA  0009  (right)  airfoils.  The  rectangular  cut-out  in  NACA  0009  airfoils  is 
used  to  mount  a  support  shaft 

The  other  airfoil  shape  considered  is  NACA  0009  with  chord  length  of  12  cm  and  aspect  ratio  of 
2.8.  A  rigid  and  flexible  versions  of  this  airfoil  are  fabricated  with  the  latter  having  the  trailing 
flexible  portion  spanning  50%  of  the  chord  length.  This  airfoil  is  used  in  experiments  in  the  full- 
scale,  61  cm  x  61  cm,  test  facility. 


IV.  3.  Velocimetry  Techniques 

The  primary  methodology  empolyed  for  measurement  of  the  velocity  field  are  Particle  Image 
Velocimetry  (PIV)  and  Molecular  Tagging  Velocimetry  (MTV).  PIV  is  used  for  measurements 
in  the  wake  of  the  long-tailed  airfoil  in  the  %-scale  facility.  These  measurements  are  done  using 
silver-coated  hollow-glass-spheres  (Potter  Industries  Inc.  -  SH400S20)  seed  with  a  mean 
diameter  of  about  20  pm.  A  500  mW  Lasiris  Magnum  SP  Laser  operating  at  a  wavelength  of 
680  nm  is  employed  to  illuminate  the  glass  spheres  in  the  wake  of  the  airfoil  at  about  mid-span. 
A  12-bit  1392  pixel  x  1024  pixel  PCO  Pixelfly  QE  camera  controlled  via  Camware  is  used  to 
record  the  PIV  snapshots.  The  camera  and  laser  pulses  are  synchronized  using  a  Stanford 
Research  Systems  DG535  digital  delay  generator.  One  thousand  pairs  of  images  are  acquired 
and  processed  for  each  case  at  the  maximal  frequency  allowed  by  our  setup  (5  Hz).  The  PIV 
region  of  interest  is  limited  to  the  wake  region  and  is  about  1.5  chord  long  in  the  streamwise  (x) 
direction  and  2  chord  wide  in  the  cross-stream  (y )  direction:  see  Figure  4.  Interrogation  of  PIV 
images  is  done  using  the  code  by  Thielicke  and  Stamhuis  (http://pivlab.blogspot.com),  yielding  a 
velocity  vector  for  every  16  pixel  x  16  pixel  interrogation  volume  (corresponding  to 
0.0283cx0. 0283c).  Vector  validation  is  implemented  using  a  standard  deviation  and  a  local 
median  check.  The  percentage  bad  vectorsis  typically  2  -  5%.  These  are  replaced  via 
interpolation. 


On  the  other  hand,  one-  (1C)  and  two-component  (2C)  MTV  measurements  are  employed  for 
velocity  measurements  in  the  0.61  m  x  0.61  m  water  tunnel.  The  lc-MTV  measurements  are 
predominantly  used  to  characterize  the  separation  bubble  around  a  stationary  SD  7003  airfoil. 
The  details  of  these  measurements  are  given  in  Olson  et  al.  (2013).  On  the  other  hand,  the  2C- 
MTV  measurements  are  utilized  in  the  experiments  focused  on  studying  the  vortex  wake 
structure  behind  rigid  and  flexible  harmonically-pitching  NACA  0009  airfoils.  For  these 
measurements,  the  beam  from  a  Coherent  COMPexPro  205  XeCl  (308  nm)  Excimer  laser  is 
directed  through  sheet  forming  optics,  split  using  a  50/50  beam  splitter  and  directed  through 
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beam  expanders  and  beam  blockers  before  passing  through  the  quartz  windows  of  the  test  section 
to  create  a  grid  pattern  with  intersection  angles  of  105°  (see  Figure  5).  Measurements  are 
performed  in  a  region  downstream  of  the  trailing  edge  and  extending  in  the  streamwise  (x) 
direction  from  0.1c  to  0.7c,  downstream  of  the  trailing  edge,  and  from  -0.2c  to  0.2c  in  the  cross¬ 
stream  (y )  direction.  The  images  are  captured  using  a  14-bit  1392  x  1024  pixel  PCO  USB 
Pixelfly  camera  paired  with  a  Nikkor  50  mm  f/1.2  lens.  The  first  image  is  acquired  -30  ps  after 
the  laser  pulse,  and  the  second  is  captured  -5.3  ms  after  the  laser  pulse.  Correlation  analysis  is 
done  on  the  two  images  using  in-house  Matlab  code  in  order  to  extract  velocity  information.  The 
resulting  time  series,  from  1400  image  pairs,  are  phase  ordered  and  then  averaged  into  64  phase 
bins  based  on  the  acquisition  (1 1.97  Hz)  and  motion  frequencies. 


PIV  measurement  plane 


2c 


1.5c 


Figure  4.  Illustration  of  the  long-tailed  airfoil  motion  trajectory  and  parameters,  and  the  PIV 
measurement  domain:  U ,  is  the  freestream  velocity,  a0  is  the  pitch  angle  oscillation  amplitude,  hp  is  the 
rigid  airfoil’s  trailing  edge  peak-to-peak  deflection  in  the  cross-stream  (y)  direction,  and  hpaclual  is  the 
corresponding  deflection  for  the  flexible  airfoil 
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Figure  5.  An  image  of  a  grid  of  tagged  fluid  molecules  in  the  wake  of  a  rigid  NACA  0009  airfoil 
undergoing  harmonic  pitch  oscillation  around  the  14-chord  point  with  reduced  frequency  k  =  13.2  and 
Strouhal  number  St  =  0.25.  The  grid  is  imaged  with  a  long  time  delay  relative  to  the  initial  excitation  of 
the  molecules  by  UV  laser  in  order  to  visualize  the  wake  vortex  pattern 

IV.  4.  Computational  Methods 

All  simulations  are  performed  employing  the  extensively  validated  high-order  FDL3DI  Navier- 
Stokes  solver  (described  in  more  detail  in  Visbal  and  Gaitonde,  1999  &  2002),  which  was 
developed  at  the  Computational  Sciences  Branch,  AFRL  Wright-Patterson  Air  Force  Base.  In 
this  code,  a  finite-difference  approach  is  used  to  discretize  the  governing  equations,  and  all 
spatial  derivatives  are  obtained  employing  a  6th-order  compact-differencing  scheme.  In  order  to 
eliminate  high-frequency  spurious  components,  an  8th-order  Pade-type  low-pass  spatial  filtering 
operator  (Visbal  and  Gaitonde,  1999  and  2002)  is  also  incorporated.  This  filter  is  applied  to  the 
conserved  variables  along  each  transformed  coordinate  direction  after  each  time  step  or  sub¬ 
iteration.  For  transitional  and  turbulent  flows,  this  filtering  technique  provides  an  effective 
implicit  LES  approach. 

Initially,  an  O-mesh  grid  is  used  for  the  computation.  However,  as  demonstrated  in  Figure  6,  the 
grid  is  found  to  be  coarse,  causing  superficially  fast  decay  of  vorticity  in  the  airfoil  wake.  To 
remedy  this  problem,  an  overset  grid  with  higher  spatial  resolution  is  used  in  the  wake,  yielding 
vorticity  decay  which  agrees  with  the  theoretical  decay  of  a  two-dimensional  vortex  having 
Gaussian  core  vorticity  distribution. 
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New  Overset  Grid  {Q-mesh:  651x392x3) 


Figure  6.  Demonstration  of  the  effect  of  using  an  overset  grid  for  computing  the  flow  in  the  wake  of  a 
NACA  0012  airfoil  undergoing  harmonic  pitch  oscillation  around  %-chord  point  with  aQ  =  2°,  k  =  6.0, 
Rec  =  12,000,  andM=  0.1:  non-dimensional  vorticity  ( coc/U^ )  maps  with  (top)  and  without  (bottom)  an 
overset  grid.  Flow  is  from  left  to  right. 


Since  the  simulation  code  is  based  on  the  governing  equations  for  compressible  flow,  attention  is 
also  paid  to  the  effect  of  the  Mach  number  (M)  on  the  simulation  results.  It  is  found  that  M  could 
have  a  significant  influence  on  the  wake  structure  even  for  M  values  of  order  0. 1  and  less.  A 
particularly  striking  example  is  shown  in  Figure  7  for  asymmetric  oscillations  of  a  NACA  0012 
airfoil.  As  seen  from  the  figure,  the  vortex  structure  in  the  far  wake  consists  of  only  two 
opposite-sign  vortices  per  wavelength  for  M  =  0.075  and  higher.  This  is  inconsistent  with 
experimental  observations  (also  depicted  in  Figure  7),  which  show  each  wavelength  to  contain  a 
positive  vortex  near  the  wake  centerline  and  two  negative  vortices  on  opposite  sides  of  the 
positive  vortex.  This  type  of  structure  could  only  be  captured  in  the  computations  for  M=  0.05  or 
less. 
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Figure  7.  Demonstration  of  the  effect  of  Mach  number  on  the  computed  wake  vortex  structure  of  the 
flow  in  the  wake  of  a  NACA  0012  airfoil  undergoing  asymmetric  pitch  oscillation  around  14-chord 
point  with  a0  =  2 °,k=  6.68  and  Rec  =  12,000:  angle  of  attack  variation  (top  right)  exhibits  faster  pitch- 
up  (38%  of  oscillation  period)  than  pitch-down  motion,  computed  non-dimensional  vorticity  ( coc/Ux ) 
maps  (left)  for  different  Mach  numbers  (with  Mach  number  increasing  from  bottom  to  top),  and  dye 
flow  visualization  (bottom  right).  Flow  is  from  left  to  right  for  experimental  and  computational  results. 

IV.  Results 


IV.  1.  Low-Reynolds-Number  Separation  over  SD  7003  Airfoil 

Coordinated  computational  and  experimental  efforts  are  undertaken  to  investigate  the  separation 
bubble  associated  with  the  low  Reynolds  number  flow  (Rec  =  2x1 04-  4><1()4)  around  an  SD7003 
airfoil  (Figure  8)  at  different  angles  of  attack.  Previous  work  has  already  demonstrated  a  great 
deal  of  variability  in  the  separation  and  reattachment  locations  (Xs  and  Xr,  respectively) 
measured  under  the  same  nominal  flow  conditions  but  in  three  different  test  facilities.  In  the 
present  work,  the  difficulty  in  accurately  determining  the  separation  and/or  reattachment 
locations  of  separation  bubbles  over  airfoils  at  low  Reynolds  number  is  further  underlined  where 
even  results  from  the  coordinated  experiments  and  computations  at  the  same  Reynolds  number 
and  angle  of  attack  can  exhibit  substantial  differences:  see  Figure  9.  By  methodically  examining 
the  experimental  uncertainty  in  determination  of  the  separation  and  reattachment  locations,  we 
are  able  to  identify  various  factors  which  render  accurate  characterization  of  the  separation 
bubble  at  low  Reynolds  number  challenging.  These  factors  are:  (I)  the  ability  to  determine  the 
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angle  of  attack  accurately,  (II)  the  sensitivity  of  the  flow  to  relatively  small  Reynolds  number 
variation  under  certain  conditions,  (III)  the  required  near-wall  measurement  spatial  resolution, 
and  (IV)  the  variation  in  freestream  turbulence  intensity  from  one  test  facility  to  another.  Details 
of  this  work  can  be  found  in  Olson  (201 1),  Olson  et  al.  (201 1)  and  Olson  et  al.  (2013). 


-a 


Figure  8.  Illustration  of  the  SD7003  airfoil  shape,  along  with  ten  overlapping  fields  of  view  used  to  capture 
the  near  surface  flow  using  high  resolution  lc-MTV  measurements.  Additional  fields  of  view  are  added 
above  these  depicted  to  capture  the  entire  shear  layer,  as  needed,  and  the  freestream  conditions 
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Figure  9.  Comparison  of  separation  and  reattachment  locations  obtained  from  different  studies  of  the  flow 
over  an  SD7003  airfoil  at  a  nominal  Rec  =  4*  104.  A  square  symbol  is  used  for  reattachment  location  whereas 
a  circle  is  used  for  separation  location  data.  Splines  are  used  to  connect  points  to  highlight  the  trends. 
Current  results  are  obtained  in  two  different  ways:  (I)  from  the  first  reliable  velocity  measurement  above  the 
surface;  (II)  using  a  polynomial  fit  to  the  reliable  data  and  the  no-slip  boundary  condition  at  the  wall. 


IV.  2.  Effect  ofTE  Flexibility  on  Wake  Vortex  Pattern 

IV.2.a.  Studies  of  the  long-tailed  airfoil  model.  In  these  studies,  PIV  measurements  are 
conducted  in  the  wake  of  harmonically  pitching  airfoils  possessing  a  flexible  ‘tail’  (see  Figure  4 
and  associated  discussion),  the  structural  flexibility  of  which  is  varied  over  three  orders  of 
magnitude.  Data  are  obtained  for  different  oscillation  amplitudes  and  frequencies  and  chord 
Reynolds  number  of  Rec  =  1950.  The  focus  of  the  data  analysis  is  on  extracting  the  circulation 
(T)  and  streamwise  and  cross-stream  spacing  (a  and  b  respectively)  of  the  vortices  shed  into  the 
wake:  see  Figure  10  for  definitions.  The  influence  of  the  airfoil  flexibility  on  how  these 
parameters  vary  with  change  in  the  amplitude  and  frequency  of  oscillation  is  investigated. 
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The  results  show  that  at  a  given  reduced  frequency,  an  optimum  airfoil  flexibility  exists  which 
results  in  the  earliest  switch  in  the  wake  vortex  pattern  from  classical  to  reverse  von  Karman 
street  (i.e.  b  switching  sign  from  negative  to  positive)  with  increasing  oscillation  amplitude.  This 
optimum  flexibility  shifts  towards  stiffer  tails  with  increasing  reduced  frequency.  More 
interestingly,  a  good  correlation  for  all  airfoils  examined  is  found  for  the  variation  of  b/c  with 
Strouhal  number  (St  =fhp/ f/x;/is  oscillation  frequency  and  hp  is  as  defined  in  Figure  4)  if  the 
latter  is  defined  using  the  actual  oscillation  amplitude  of  the  trailing  edge  (Stte  =  fhp,actuai  /  £/«,; 
hp, actual  is  as  defined  in  Figure  4).  These  results  are  depicted  in  Figure  11.  Similar  correlations  are 
also  found  for  the  dependence  of  the  circulation,  streamwise  spacing  and  convection  velocity 
when  using  the  latter  definition  of  the  Strouhal  number.  Details  may  be  found  in  Monnier  et  al. 
(2013). 


Figure  10.  Definition  of  geometrical  parameters  of  the  vortex  pattern  in  the  wake  of  the  long-tailed  airfoil 
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Figure  11.  Dependence  of  b/c  on  St  (left)  and  Stte  (right)  for  the  long-tailed  airfoil.  Colors  represent  airfoil 
flexibility:  blue:  rigid,  red:  G\  =  5.6,  green:  G\  =  1.4,  cyan:  G\  =  0.7.  Symbols  represent  amplitude  aQ\ 
circle:  2°,  right  triangle:  3°,  square:  4°,  star:  5°,  left  triangle:  6°,  diamond:  7°.  Symbol  filling  represents 
frequency:  open  symbols:  k  =  4.89,  filled  symbols:  k  =  6.19,  filled  symbol  with  black  outline:  k  =  7.74.  b/c 
uncertainty  is  estimated  to  be  ±0.028c  (corresponding  to  ±  one  PIV  grid  spacing). 


IV.2.b.  Preliminary  studies  of  the  NACA-0009  airfoil  model.  This  component  of  the 
investigation  is  still  ongoing.  The  results  shown  here  represent  some  preliminary  observations  of 
differences  between  the  wake  vortex  pattern  observed  in  the  wake  of  a  rigid  and  flexible  NACA- 
0009  airfoils,  obtained  using  2C-MTV  measurements  at  Rec  =  12,500.  For  the  case  examined 
both  airfoils  undergo  pitch  oscillation  around  %-chord  point  with  an  amplitude  of  2.25  degrees 
and  reduced  frequency  of  k  =  13.2  (St  =  0.25).  The  trailing  edge  (TE)  displacement  of  both 
airfoils  is  measured  from  the  MTV  images  and  the  results  are  shown  in  Figure  12.  The  flexible 
airfoil  TE  movement  is  found  to  have  27%  larger  amplitude  than  the  rigid  airfoil  and  to  exhibit  a 
phase  lag  of  34  degrees.  The  vorticity  field  for  both  airfoils  is  shown  in  Figure  13  at  the  same 
phase  of  the  oscillation  cycle.  Using  these  results,  parameters  of  the  wake  vortex  pattern  are 
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computed  and  compared  in  Table  1.  As  can  be  seen  from  the  table,  the  airfoil  flexibility  results  in 
larger  cross-stream  vortex  spacing  and  vortex  circulation  (even  though  the  imposed  airfoil 
motion  in  both  cases  is  the  same).  These  differences  (i.e.  increased  b/c  and  T)  ,  which  are  clearly 
related  to  the  different  boundary  condition  associated  with  the  deformation  of  the  flexible  airfoil 
(e.g.  the  difference  in  TE  motion  seen  in  Figure  12),  are  expected  to  result  in  producing  larger 
mean  thrust  on  the  airfoil.  In  fact,  the  mean  streamwise  force  coefficient  (Cf)  is  computed  for 
these  two  cases  using  the  momentum  integral  equation  (see  Bohl  and  Koochesfahani  (2009)  for 
the  details),  and  the  results  are  shown  in  Figure  14  along  with  data  from  an  earlier  study  by  Bohl 
(2002),  which  covered  a  substantially  wider  range  of  oscillation  frequencies.  As  can  be  seen 
from  the  figure,  the  flexible  airfoil  does  in  fact  results  in  a  larger  mean  thrust  coefficient 
(specifically  0.105  in  comparison  to  0.046:  a  2.3  fold  increase).  It  is  also  interesting  to  note  how 
the  current  Cf  data  along  with  data  from  Bohl  (2002)  seem  to  collapse  on  a  single  trend  line 
when  plotted  versus  the  Strouhal  number  based  on  the  actual  TE  oscillation  amplitude.  These 
results  are  consistent  with  the  conclusions  made  in  section  IV.2.a.  regarding  the  dependence  of 
the  wake  vortex  parameters  of  the  long-tailed  airfoil  with  different  tail  flexibilities  on  the 
Strouhal  number. 


Figure  12.  Measured  trailing  edge  deflection  of  the  rigid  (red  trace)  and  flexible  (green  trace)  NACA-0009 
airfoils  over  one  oscillation  cycle.  Deflection  is  shown  in  units  of  “milli-chord”  (me) 


Figure  13.  Velocity  field  superposed  on  top  of  the  vorticity  field  measured  in  the  wake  of  the  rigid 
(left)  and  flexible  (right)  NACA-0009  airfoils  at  the  same  phase  of  the  oscillation  cycle.  The  color 
bar  indicates  the  vorticity  magnitude  in  s'1,  and  the  arrow  above  the  top-left  corner  of  each  plot 
represents  the  freestream  velocity  of  20  cm/s 
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Table  1.  Comparison  between  the  wake  vortex  parameters  for  the  rigid  and  flexible  NACA-0009  airfoils 


Parameter 

Rigid  Airfoil 

Flexible  Airfoil 

Difference 

b/c 

0.068 

0.088 

28% 

Peak  Vorticity  (s'1) 

75 

91 

21% 

f(cm2/s) 

52 

68 

30% 

Stg 

Figure  14.  Comparison  of  the  mean  streamwise  force  coefficient  obtained  from  the  current  results 
(Olson)  with  data  obtained  by  Bohl  (2002)  in  the  wake  of  NACA-0012  airfoil  exhibiting  harmonic 
pitch  oscillation  around  the  %-chord  point 


IV.  3.  Reynolds  Number  Effects 

Our  low  Reynolds  number  data  at  Rec  «  2,000  (Monnier,  et  al.  2013)  indicate  that  the  crossover 
reduced  frequency  for  the  vortex  arrangement  switch  (i.e.  from  wake  to  jet  profile)  occurs  at  k  « 
7.8  for  rigid  airfoil  at  amplitude  of  A  =  2°.  This  is  significantly  higher  than  k  ~  5.9  previously 
measured  as  the  crossover  condition  at  Rec  ~  12,600  for  the  same  oscillation  amplitude  (Bohl  and 
Koochesfahani,  2009).  As  we  will  show  below,  similar  Reynolds  number  effects  also  occur  in 
connection  to  the  conditions  needed  for  thrust  generation. 

To  examine  the  Reynolds  number  effects,  CFD  studies  were  carried  out  using  the  FDL3DI  solver 
of  Visbal  at  the  AFRL.  Figure  15  depicts  the  computed  crossover  reduced  frequency  k  for  switch 
in  vortex  arrangement  in  comparison  with  existing  limited  data  from  both  experiments  and  other 
computations.  A  clear  Reynolds  number  dependence  is  revealed  in  this  figure.  We  note  the  very 
good  agreement  between  the  crossover  k  from  our  current  CFD  results  and  our  previous 
experiments.  However,  notable  differences  from  previous  computional  results  are  also  observed. 
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Figure  15.  Reynolds  number  dependence  of  crossover  reduced  frequency  k  for  the  vortex 
arrangement  switch  (NACA  0012;  Oscillation  amplitude  T  =  2°). 


The  computed  influence  of  Reynolds  number  on  the  crossover  reduced  frequency  k  for  the  drag- 
to-thrust  switch  is  shown  in  Figure  16.  Consistent  with  the  previous  finding  of  Bohl  & 
Koochesfahani  (2009),  the  values  of  k  for  the  drag-to-thrust  switch  are  higher  than  those  for 
switch  in  vortex  arrangement.  Note  that  the  drag-to-thrust  crossover  reduced  frequency  decreases 
by  about  a  factor  of  two  (from  13.91  to  6.86)  as  the  Reynolds  number  increases  from  2,000  to 
22,000.  For  reference,  we  would  mention  that  the  Euler  calculations  of  Ramamurti  &  Sandberg 
(2001)  estimated  the  drag-to-thrust  crossover  to  be  at  k  «  3.35  in  the  inviscid  limit,  significantly 
lower  than  the  range  of  values  for  the  viscous  cases  studied  here. 

Figure  17  illustrates  the  comparison  of  our  data  with  those  in  the  literature  in  terms  of  the 
crossover  Strouhal  number  St  for  the  drag-to-thrust  switch  versus  Reynolds  number.  Using  the 
Strouhal  number,  which  is  based  on  the  peak-to-peak  trailing  edge  amplitude,  allows  us  to 
compare  results  across  different  pitch  oscillation  amplitudes.  The  results  in  Figure  17  are 
instructive  in  two  ways.  First,  our  computed  results  show  an  excellent  quantitative  agreement 
with  our  previous  experiment  (Bohl  and  Koochesfahani,  2009)  and  also  the  estimate  provided  by 
our  vortex  array  model  (Naguib,  et  al.,  2011)  applied  to  incomplete  data  of  Koochesfahani 
(1989).  Second,  there  is  significant  disagreement  among  the  various  computational  studies  in 
their  prediction  of  the  crossover  frequency,  in  some  cases  by  nearly  a  factor  of  two. 
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Figure  16.  Evidence  for  Reynolds  number  dependence  of  crossover  reduced  frequency  for  the 
drag-to-thrust  switch  (NACA  0012;  Oscillation  amplitude  A  =  2°). 


We  note  that  the  quantitative  agreement  in  Figure  17  between  our  current  CFD  results  and  data 
from  our  previous  experiments  required  significant  attention  to  numerical  issues  (resolution,  grid 
generation,  convergence,  etc.).  A  major  difficulty  in  assessing  the  adequacy  of  the  various 
computational  studies  of  unsteady  low-Re  airfoils  have  been  a  lack  of  detailed  quantitative  data 
from  experiments  with  conditions  matched  to  the  computations,  and  comparisons/validations  are 
often  made  on  a  qualitative  basis.  Qualitative  comparisons  seem  to  be  prevalent  even  when 
appropriate  experimental  data  are,  in  fact,  available. 

The  Reynolds  number  effects  shown  in  Figures  15-17  have  not  been  previously  addressed,  to  our 
knowledge,  but  they  have  important  implications.  They  may  play  a  partial  role  in  the 
disagreements  that  are  routinely  noted  between  CFD  results  and  experiments.  More  importantly, 
they  indicate  that  the  airfoil  operating  condition  needed  for  generation  of  thrust  is  unique  to  the 
particular  value  of  chord  Reynolds  number. 

The  results  of  this  effort  are  being  packaged  in  a  manuscript  for  journal  submission. 

IV.  4.  Development  ofMTV-based  Wall-Shear  and  -Pressure  Measurement  Technique 

Evaluation  of  the  Navier-Stokes  equations  at  the  surface  of  a  moving  solid  object  leads  to  the 
following  expression  for  the  wall  pressure  distribution  along  the  surface: 
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Figure  17.  Evidence  for  Reynolds  number  dependence  of  crossover  Strouhal  number  for  the 
drag-to-thrust  switch  (NACA  0012). 
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In  this  expression  r/,  £  are  coordinates  normal  and  tangent  to  the  surface,  pw(  g0)  is  a  “boundary 
condition”  specifying  the  pressure  at  an  arbitrarily  selected  point  on  the  wall,  aw.^  is  the  surface- 


tangent  component  of  the  wall  acceleration,  u  ,  is  the  flow  velocity  component  in  the  £  direction 

and  p  is  the  dynamic  viscosity.  Inspection  of  Equation  [2]  shows  that,  given  a  prescribed 
airfoil/wing  motion  trajectory  and  measurement  of  the  wall  pressure  at  an  arbitrary  single 
reference  point  on  the  surface  (e.g.  at  the  leading  edge),  the  pressure  distribution  around  the 

d2ut 

airfoil,  at  a  given  spanwise  location,  can  be  obtained  from  the  measurement  of  2  over  the 


dr/ 


circumference  of  the  airfoil/wing. 


In  our  work,  we  take  advantage  of  the  high  spatial  (pixel  level)  resolution  of  the  single¬ 
component,  line-tagging  mode  of  MTV  (lc-MTV)  for  near- wall  measurements.  The  near- wall 
capability  is  due  to  both  the  molecular  nature  of  MTV  and  also  to  little/no  wall  glare.  By 
orienting  the  line  tags  normal  to  the  surface  (a  condition  imposed  for  convenience  but  can  be 
relaxed),  lc-MTV  essentially  provides  a  “continuous”  distribution  of  in  the  rj  direction.  Thus, 
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(after  correction  for  the  wall-normal 


the  curvature  of  the  line  tags  at  the  wall  is  physically 


d2uf 


dr] 2 


velocity  due  to  the  motion  of  the  airfoil  as  mentioned  above).  In  addition  to  measurements  of  the 
surface  pressure,  the  line-tagging  MTV  mode  can  be  used  to  obtain  the  £  component  of  the  wall 
shear  stress  tw-^.  In  this  case,  the  slope  of  the  line  tags  at  the  wall  will  be  used  to  calculate  the 
shear  stress  using: 
du, 

P] 


dr] 


An  initial  study  was  carried  out  to  assess  the  accuracy  of  computing  the  second  derivative  of  the 
velocity  profile  at  the  wall  from  MTV  data  and  how  it  depends  on  the  spatial  resolution  of  the 
measurement,  the  streamwise  pressure  gradient,  the  boundary  layer  thickness,  and  the  order  of 
accuracy  of  the  numerical  derivative  scheme.  This  study  utilized  the  Falkner-Skan  boundary 
layer  solutions  as  a  model  for  boundary  layers  developing  under  non-zero  pressure  gradient, 
along  with  the  CFD  solution  of  the  flow  over  an  airfoil  (SD7003)  for  an  analysis  of  a  complete 
flow  field.  The  results  of  this  work  confirmed  the  feasibility  of  surface  pressure  estimation  from 
multi-line  MTV  measurements  and  also  established  the  various  constraints  that  impact  the 
accuracy  of  the  approach. 

While  the  line-tagging  mode  of  MTV  provides  highly  resolved  measurements  of  a  single 
component  of  the  velocity  normal  to  the  line  tags,  the  accuracy  of  the  measured  velocity  can 
suffer  from  the  “V-bias”  problem.  This  problem  is  connected  with  the  component  of  the  velocity 
that  might  exist  along  the  line  tags,  causing  an  ambiguity  in  uniquely  determining  the 
displacements  of  the  ‘individual  fluid  particles’  comprising  the  ( continuous )  tagged  line.  A 
detailed  analysis  of  the  error  led  to  a  novel  correction  method  involving  the  acquisition  of  MTV 
images  at  multiple  time  delays.  The  effectiveness  of  this  ‘multi-time-delay’  approach  for 
correcting  the  lc-MTV  data  has  been  validated  experimentally  and  the  results  appear  in  a  recent 
journal  publication  (Hammer,  et  al.,  2013) 

We  have  recently  carried  out  a  quantitative  experimental  demonstration  of  the  MTV-based 
capability  for  measuring  surface  pressure  and  shear  stress  distributions.  The  experiments 
considered  the  flow  around  a  stationary  circular  cylinder  in  a  water  tunnel  at  a  cylinder  Reynolds 
number  (based  on  diameter)  of  Re  =  6,000.  For  this  flow  field,  Equation  [2]  reduces  to  the 
following  expression. 


Pw(%)-Pw(£o) 


\d£, 


[4] 


Figure  18  shows  a  composite  image  of  the  tagged  regions  acquired  7  ms  after  initial  tagging.  The 
laser  beam  was  arranged  to  exit  the  surface  of  the  transparent  circular  cylinder  in  a  radial 
direction  at  a  series  of  angular  locations.  The  wall  normal  spatial  resolution  was  37  pm,  resulting 
in  about  27  points/mm  in  the  velocity  profile.  Two  samples  of  the  measured  mean  velocity 
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profiles  are  indicated  in  this  figure  for  the  two  selected  angular  positions.  The  highly  resolved 
nature  of  the  velocity  data  near  the  surface  is  noted. 


Figure  18.  Composite  delayed  image  of  the  radially  tagged  regions  around  the  circular  cylinder  at  Re  =  6,000. 
The  delayed  image  is  acquired  7  ms  after  initial  tagging.  Wall  normal  resolution  is  37  pm  (27  points/mm). 
The  two  insets  show  the  measured  velocity  profiles  at  the  two  particular  angular  locations  marked. 


A  3rd  order  polynomial  fit  to  near  wall  data  was  used  to  estimate  the  second  and  first  order 
derivatives  of  the  velocity  at  the  surface.  The  second  order  derivative  was  used  in  conjunction 
with  Equation  [4]  to  determine  the  surface  pressure  distribution,  along  with  specifying  the 
boundary  condition  at  the  front  stagnation  point.  Results  are  shown  in  Figure  19  in  terms  of  the 
distribution  of  surface  pressure  coefficient  Cp.  We  can  see  that  our  results  at  Re  =  6,000  are  in 
very  good  agreement  with  published  results  at  neighboring  values  of  Reynolds  number. 

Comparison  of  the  measured  surface  shear  stress  distribution  with  literature  data  is  more 
challenging  since  this  variable  is  often  not  included  in  publications.  However,  we  were  able  to 
obtain  from  the  Visbal  group  at  AFRL  the  shear  stress  data  for  the  computations  reported  in 
Rizzetta  &  Visbal  (2009)  for  Re  =  10,000.  To  scale  their  shear  stress  data  to  the  Reynolds 
number  of  our  experiment  (Re  =  6,000),  we  note  that  for  zero  pressure  gradient  (Blasius) 
boundary  layers  the  surface  friction  coefficient  C/  scales  as  1/ VRe  .  This  scaling  also  holds  for 
non-zero  pressure  gradient  (Falkner-Skan)  boundary  layer  solutions.  Since  according  to  Figure 
19  the  pressure  distribution  is  very  weakly  dependent  on  the  Reynolds  number,  over  the  range  Re 
=  6,000  -  10,000,  the  boundary  layer  evolves  under  very  similar  external  pressure  gradient  and 
we  expect  the  1/  VRe  scaling  for  C/  to  continue  to  hold.  The  measured  surface  friction  coefficient 
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Cf  and  that  arrived  at  from  the  computations,  using  the  scaling  just  described,  are  shown  in 
Figure  20.  We  point  out  the  excellent  agreement  between  the  data  sets  in  this  figure. 

The  results  of  this  development  are  being  presented  at  the  2013  meeting  of  the  APS/DFD  and  a 
journal  manuscript  is  also  being  prepared. 


Figure  19.  Comparison  of  the  measured  circular  cylinder  surface  pressure  coefficient  against  literature  data. 


Figure  20.  Comparison  of  the  measured  circular  cylinder  surface  shear  stress  distribution  against  the 
computational  results.  Scaling  according  to  1/  Vlte  was  used  to  transform  the  computational  results  at 
Re  =  10,000  to  the  Reynolds  number  of  the  experiment. 
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V.  Summary 

A  coordinated  experimental  and  computational  investigation  is  carried  out  to  study  the 
fundamental  aerodynamics  of  oscillating  airfoils  in  the  chord  Reynolds  number  range  of  2x1 03  - 
6xl04.  The  objective  of  the  investigation  is  to  establish  the  connection  between  the  airfoil  motion 
trajectory,  the  trailing  edge  flexure,  the  time  history  of  vorticity  flux  at  the  trailing  edge,  the 
pattern  of  shed  vorticity  and  its  evolution,  and  the  conditions  for  generation  of  thrust/drag  and 
lift.  The  following  list  summarizes  the  key  findings  of  the  study: 

•  In  computing  the  flow  over  unsteady  airfoils,  special  attention  must  be  paid  to  the 
computational  grid  resolution  if  one  is  to  study  the  characteristics  of  the  wake  vortex  pattern. 
The  resolution  requirement  in  the  wake  could  be  more  stringent  than  required  to  capture  the 
force  on  the  airfoil  without  worrying  about  the  connection  to  the  wake  behavior.  In  the 
present  study,  the  grid  resolution  requirements  in  the  wake  are  met  using  a  fine  Cartesian  grid 
superposed  on  top  of  an  O  grid  in  the  wake  of  the  airfoil. 

•  When  using  a  compressible-flow  code  to  compute  the  flow  over  unsteady  airfoils, 
compressibility  may  have  substantial  influence  on  the  results  even  at  Mach  numbers  as  low 
as  0. 1  or  lower.  The  influence  is  particularly  significant  at  high  reduced  frequency.  A  Mach 
number-independent  solution  as  M— >0  (or  one  that  is  consistent  with  experimental 
observations)  should  be  sought  systematically  if  the  target  application  is  that  involving 
incompressible  flow. 

•  Experimental  determination  of  the  characteristics  of  separation  bubbles  over  airfoils  at  low 
Reynolds  number  is  challenging.  The  factors  contributing  to  this  difficulty  have  been 
identified  and  examined  systematically  in  the  present  work. 

•  The  wake  vortex  characteristics  and  the  mean  thrust  coefficient  of  rigid  and  flexible  airfoils 
are  similar  if  characterized  using  the  Strouhal  number  based  on  the  actual  oscillation 
amplitude  of  the  trailing  edge,  rather  than  using  the  Strouhal  number  based  on  the  imposed 
motion. 

•  For  a  given  oscillation  reduced  frequency,  an  optimum  trailing  edge  flexibility  exists,  which 
will  result  in  the  switch  in  the  vortex  wake  pattern  to  the  “reversed”  von  Karman  vortex 
street  at  the  smallest  possible  oscillation  amplitude.  This  optimum  flexibility  is  such  that  the 
higher  the  reduced  frequency  of  oscillation,  the  stiffer  the  structure  of  the  optimum  trailing 
edge  material. 

•  The  Reynolds  number  has  a  very  important,  yet  unrecognized  hitherto,  influence  on  the 
reduced  frequency  and  Strouhal  number  at  which  the  wake  vortex  pattern  assumes  the 
“reversed”  von  Karman  vortex  street  configuration  and  the  net  force  acting  on  the  airfoil 
switches  from  drag  to  thrust.  This  effect  is  documented  systematically  in  the  present  work 
and  the  underlying  physics  leading  to  this  influence  is  currently  under  examination. 

•  New  multi-time-delay-based  method  is  introduced  for  improvement  of  the  accuracy  of  high- 
resolution  (pixel  level)  one-component  Molecular  Tagging  Velocimetry  (lc-MTY).  The 
interest  in  lc-MTV  is  motivated  by  the  suitability  of  the  method  for  conducting  near- wall 
boundary-layer-resolved  measurements.  This  suitability  is  capitalized  upon  in  the  present 
investigation  in  order  to  extract  the  wall-shear-stress  and  wall-pressure  information.  The 
approach  is  validated  by  conducting  measurements  over  a  cylinder  in  cross  flow,  and 
comparing  the  results  against  existing  experimental  and  computational  studies  of  the  same 
flow.  This  approach  for  measuring  the  wall  shear  stress  and  pressure  provides  a  clever 
solution  for  measuring  these  quantities  over  flexible  and  moving  surfaces  where  such 
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measurements  are  challenging  to  obtain  using  surface  mounted  sensors,  pressure  taps,  or 
other  traditional  means.  A  particularly  significant  area  of  application  of  the  new  technique  is 
in  aeroelastic  problems  where  access  to  velocity  field  information  and  the  concurrent  stresses 
acting  on  a  moving  structure  would  yield  tremendous  insight. 
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